Critical behavior of the Coulomb-glass model in the zero-disorder limit: 
Ising universality in a system with long-range interactions 
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The ordering of charges on half-fiUed hypercubic lattices is investigated numericaUy, where elec- 
troneutrality is ensured by background charges. This system is equivalent to the s = 1/2 Ising lattice 
model with antiferromagnetic 1/r interaction. The temperature dependences of specific heat, mean 
staggered occupation, and of a generalized susceptibility indicate continuous order-disorder phase 
transitions at finite temperatures in two- and three-dimensional systems. In contrast, the suscepti- 
bility of the one-dimensional system exhibits singular behavior at vanishing temperature. For the 
two- and three-dimensional cases, the critical exponents are obtained by means of a finite-size scaling 
analysis. Their values are consistent with those of the Ising model with short-range interaction, and 
they imply that the studied model cannot belong to any other known universality class. Samples of 
up to 1400, 112^, and 22^ sites are considered for dimensions 1 to 3, respectively. 

PACS numbers; 64.60.F-, 05.70.Jk, 71.10.-w, 02.70.Uu 



I. INTRODUCTION 

The long range of electrostatic correlations un- 
derlies important physical effects such as screening, 
charge renormalization, charge orders and instabilities 
of plasmasi In Coulomb glasses, disordered systems of 
localized charges, electrostatic correlations induce the 
Coulomb gap in the single-particle energy spectrum:- 
The question under which conditions these glasses ex- 
hibit genuine glass transitions has been under controver- 
sial debate for a long iiuieM^^^^L^^^^^'^^ 

Lattices partially occupied by charged particles have 
been studied as models for ionic fluids, where several 
phases were identified and an Ising-type liquid-gas crit- 
ical point was observed i-'^^i-'^^'^'^i^^'^^i"'^^ In particular, by 
means of numerical simulation and finite-size scaling, 
Luijten et al. could clearly rule out this point to belong 
to any of several alternative universality classes. For a 
recent review on this field see Ref. llTi. In case of full 
occupancy, such a lattice corresponds to an ionic crystal 
and provides the zero-disorder limit of the Coulomb-glass 
problem. 

Although the ordering of ionic crystals with decreas- 
ing temperature has been studied by several groups, im- 
portant questions are still open. Analytica l^^'^^i^^ and 
simulation results^ii^ii^i^ suggest that, for the case of 
pure Coulomb interaction, staggered (antiferromagnetic) 
ordering in three-dimensional simple-cubic lattices starts 
with a continuous phase transition. According to the 
hierarchical reference theory study by Brognara et alr^ 
and the renormalization-group investigation by Ciach^i^ 
this transition should belong to the short-range Ising uni- 
versality class in spite of the long-range interaction. This 
hypothesis is suggested by the transition being the end 
of a line of Ising transitions in the lattice restricted prim- 
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itive model, which is reached for complete filling. 

However, we remind that, in the case of Ising mod- 
els with ferromagnetic long-range interactions, the phase 
transition depends on dimensionality and decay expo- 
nent of the interactions. The critical indices can vary 
from mean-field behavior to that of the short-range Ising 
universality class via intermediate valuesi ^^'^^ 

The case of antiferromagnetic long-range interactions 
is more subtle due to their inherent frustration. From the 
corresponding field-theoretic approach for a staggered 
order, one would expect that the long-distance behav- 
ior of the interactions should renormalize to short-range 
couplings However, this argument assumes validity of a 
perturbative expansion around the ground-state. Thus, 
one discounts part of the couplings between charge or 
higher-order multipole fluctuations that are associated 
with the defects of the ordered state. In the case of 
Coulomb-glass systems, the interplay of such fluctuations 
and quenched disorder is believed to underly glassy prop- 
erties. Recently, three different analytical approaches for 
the screening problem in such disordered electronic sys- 
tems have been developed, where two of them contain the 
non-disordered case as natural limit. ^'^'^'^ In this context, 
precise numerical experiments on the phase transition of 
lattice models without disorder are called for. 

The corroboration of the analytical results on the or- 
dering of ionic crystals by numerical studies is very diffi- 
cult due to the long-range interaction: Studying samples 
of up to 12^ sites, Almarza and Enciso observed their 
data for simple-cubic lattices to be consistent with the 
assumption of short-range Ising universalityj^^ However, 
they could not determine the critical exponents directly 
because the samples were still too small. A considerable 
progress was reached in Refs. ilQi and i20|, where samples 
of up to IS'^ sites were considered, so that critical expo- 
nents could be obtained by finite-size scaling. The expo- 
nent values were found to be very close to the data for 
the short-range Ising model, but this work was presented 
only in a short form. Moreover, further enlarging of the 
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sample size and diminishing of the error bars seemed de- 
sirable. Finally, Overlin et al. studied the influence of 
positional disorder at samples of up to 8^ sites and de- 
termined the critical exponent of the localization length/ 
For the limiting case of vanishing disorder, their numeri- 
cal result was consistent with both the short-range Ising 
and mean- field universality classes. 

The effects of unscreened interactions and frustration 
may become more prominent in dimensions d — 1 and 
2. However, only a few numerical studies have been de- 
voted to the behavior of such systems: In a first attempt, 
Di'az-Sanchez et al. studied samples of up to 26 and 14^ 
sites for d = 1 and 2, respectively, by means of a spin- 
glass approachjS They concluded that a phase transition 
at finite T exists only for d — 2, but not for d = 1. The 
behavior of considerably larger samples was simulated in 
Refs. and [l^, where a continuous phase transition at 
finite T was obtained for d = 2, but not for d = 1. In 
the former case, according to the obtained critical expo- 
nents, short-range Ising universality seems likely. More- 
over, Luo et al. considered two-dimensional systems with 
logarithmic interaction, corresponding to the interaction 
of homogeneously charged lines. ^'^ They obtained values 
of the critical exponents of correlation length and order 
parameter which clearly differ from the short-range Ising 
values. 

The aim of the present work is twofold: On the one 
hand, for d = 3, we numerically investigate the ordering 
in simple-cubic systems to check the analytical theories 
mentioned above. The critical exponents of correlation 
length, order parameter, specific heat, and generalized 
susceptibility are obtained by finite-size scaling. On the 
other hand, we extend these studies to the cases d = 1 
and 2. For a preliminary and less detailed version of this 
work, which was still restricted to the investigation of 
smaller systems, see Ref. [l^. 

This paper is organized as follows. In Sec. II, we intro- 
duce the model Hamiltonian including the used boundary 
conditions, and we comment on the applied numerical 
procedures. Section III presents a qualitative overview 
of the simulation results, whereas Sec. IV is devoted to 
a quantitative evaluation of the data sets by means of 
finite-size scaling. Finally, in Sec. V, the obtained crit- 
ical exponent values are discussed and compared with 
previous work. 



II. MODEL AND NUMERICAL APPROACH 

We numerically investigate the nature of the order- 
disorder transition in a system of charges considering 
the minimal model: simple-hypercubic lattices are half- 
filled with particles interacting via the Coulomb potential 
where background charges -1/2 are attached at each site 
for neutrality, 

^=^E/^^("^-l/2)K-l/2) (1) 



with 

4- = l/|r,-r,| . (2) 

Here Ua G {0, 1} denote the occupation numbers of states 
localized at sites within a d-dimensional hypercube of 
size L'^. Elementary charge, lattice spacing, dielectric 
and Boltzmann constants are all taken to be 1, so that 
the temperature T is dimensionless. Due to particle-hole 
symmetry, canonical and grand canonical Hamiltonians 
equal each other in this case. When substituting Si for 
Ui — 1/2, Eqs. (1) and (2) take the form of the s = 1/2 
Ising lattice model with Coulomb interaction of antifer- 
romagnetic character. 

For reducing finite-size effects, we impose periodic 
boundary conditions for d = I and 2, modifying fij ac- 
cording to the minimum image convention.'^''' For c? = 3, 
the same approach would give rise to an unphysical fea- 
ture: The groundstate would be a layered arrangement 
of charges instead of the expected NaCl structure in case 
L is a multiple of fouri^ Therefore, similar to Ref. [26l . 
we consider the sample to be surrounded by 26 equally 
occupied cubes in this case^ 

^ 26 ^ 1 

= + £^|r.-r,-R..| - ■ 

Here, Rj, denotes the shift of the neighboring cube k 
with respect to the central cube. Similar to the mini- 
mum image convention, the correction terms in Eq. (3) 
efficiently reduce the largest finite-size effect, namely the 
difference between the surroundings of sites close to the 
center of the cube and of sites close to its surface. Com- 
pared to the implementation of periodic boundary condi- 
tions combined with an Ewald summation, our method 
has the advantage that it does not introduce the artificial 
long-range correlations arising from the series of periodic 
images of the cell. However, in the limit L — > cxd, both 
approaches should yield the same results, see the com- 
parison with Ref. 7 in Sec. IV. 

To obtain ensemble averages of various observables, we 
follow the Metropolis approach and substitute temporal 
for ensemble averaging. Since only equilibrium proper- 
ties are of interest here, we are free to choose the dynam- 
ics so that the simulation effort is minimized. A cluster 
Monte Carlo algorithm seems not to be available for the 
antifcrromagnetic Coulomb interaction because of frus- 
tration. Thus we select the system modifications to be 
taken into account by hand: We include one-electron ex- 
change with the "surroundings" , one-electron hops over 
distances below a certain bound, and two-electron hops 
changing the occupation of four neighboring sites. The 
maximum permitted distance of one-electron hops is en- 
larged when T is lowered. 

At high T, we use the original Metropolis methodi^ 
But at low T, we take advantage of the hybrid proce- 
dure proposed in Ref. [H, which much accelerates the 
computations: Similar to the n-fold way algorithm^SS we 
deterministically calculate the rates of the transitions to 
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FIG. 1: (Color online) Temperature dependences of the spe- 
cific heat, c(r), for dimensions d = 1 to 3 as obtained from 
simulations of samples of L'' sites, d = 1: L = 100 (+), 280 
(o), and 700 (x); d = 2: L = 20 (+), 34 (o), and 58 (x); 
d = 3: L = 8 (+), 12 (o), and 18 (x). Only a part of the 
data points forming the curves is marked by symbols, as well 
as in Figs. 2 to 4, 6, and 8. The error bars are considerably 
smaller than the symbol size. 



all multiparticle states, which are accessible from the cur- 
rent state by means of a single system modification. Thus 
the dwelling time at the current state can be determined 
directly, and only one Monte Carlo step is needed per sys- 
tem modification. Moreover, the hybrid procedure con- 
nects the deterministic evaluation of weighted sums over 
all states within a low-energy subset of the configuration 
space with Monte Carlo sampling of the complementary 
high-energy subset. These two ideas had proved to be 
very efficient in studying the specific heat of the Coulomb 
glass at low temperaturesi^ 

For an efficient error control, we decompose the sim- 
ulation time considering 100 intervals with integration 
time r instead of one interval of length 100 r. In detail, 
we increase r step by step performing 50 runs for every 
T value. In each of these runs, after starting randomly 
from one multiparticle state out of a set of previously tab- 
ulated low-energy states, the sample is first equilibrated 
during a time interval t/3. Then the evolution over two 
successive time intervals r is emulated. The obtained en- 



FIG. 2: (Color online) Temperature dependence of the aver- 
age absolute value of the mean staggered occupation, {\<t\){T), 
which is defined by Eq. (5) for d = 1 to 3. For the meaning 
of the symbols o, and x see caption of Fig. 1; • marks the 
extrapolation L ^ oo explained in the text. 



sembles of 100 measuring data for each observable are 
used to estimate mean values and their statistical un- 
certainties, and to check equilibration. Based on these 
results, it is decided, whether the iteration process can 
be stopped, or whether t has to be increased further. 



III. QUALITATIVE RESULTS 

We now turn to the qualitative behavior of specific 
heat, order parameter, and susceptibility. The specific 
heat c was obtained from energy fluctuations utilizing 



(4) 



see e.g. Ref. ISOl Figure 1 shows its T and L dependences: 
For d — 2 and 3, sharp peaks of increasing height evolve 
within a small T region as L grows. Away from the peaks, 
within the T intervals presented, c is almost independent 
of L. However, for d = 1, there are only broad rounded 
peaks with L-independent height - a logarithmic T scale 
is used for d 1, in contrast to the linear scales for 
d = 2 and 3, which display far smaller T intervals. For 
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FIG. 3: (Color online) Temperature dependence of the sus- 
ceptibility, x{T), which is related to the mean staggered oc- 
cupation a by Eq. (6), for d = 1 to 3. For the meaning of the 
symbols see caption of Fig. 1. 
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FIG. 4; (Color online) Temperature dependence of the Binder 
parameter, Q{T), which is related to the mean staggered oc- 
cupation a by Eq. (7), for d = 1 to 3. For the meaning of the 
symbols see caption of Fig. 1. 



d = 1, finite-size effects are restricted to low T where the 
reliability bound decreases as L grows, compare Ref. [2^. 

Analogously to an antifcrromagnet, the order inher- 
ent in a charge arrangement n.i can be characterized 
by the mean staggered occupation a relating to a NaCl 
structured For d = 3, 

'^=ilE(-l)''^''^'"(2n,-l) (5) 

i 

where Xi, yi, and Zi denote the (integer) components 
of r^. Thus we consider the ensemble average of the 
absolute value of the mean staggered occupation {\a\) 
as order parameter. T and L dependences of (|cr|) are 
shown in Fig. 2. For d = 1, a rapid decrease of {\a\) 
with increasing T occurs already clearly below the tem- 
perature of maximum c (same T scales in Figs. 1 and 
2). This marked decrease shifts to lower T with in- 
creasing L. For d — 2 and 3, a qualitatively different 
behavior is found: {\a\) decreases rapidly just in that 
T region where the peak of c(T) evolves, and the T 
interval of rapidly diminishing {\a\) shrinks as L rises. 
In these cases, the extrapolation L ^ oo by means of 
(|cr|)(T,L) = {\a\){T, oo) + A{T)/L'^/^ seems natural 



Based on the {\a\){T,L) data for L = 34 and 58 in case 
d = 2, and for L = 12 and 18 in case d = 3, it yields 
almost sharp transitions. However, this extrapolation is 
of limited accuracy in the immediate vicinity of the tran- 
sition. 

The generalized susceptibility x, which is related to 
the response of {\a\) to a staggered field, is given by 

x-L'^((a2)-(H)2)/r, (6) 

see Ref. HH. Figure 3 shows T and L dependences of 
X- On the one hand, for d = 1, a broad peak of x(T) 
evolves with increasing L where Tmax, the temperature 
of maximum decreases. On the other hand, for d = 2 
and 3, as L rises, a narrow peak grows in just that T 
region where c(T, L) has such a feature. 

Hence, according to the described behavior of c(T, L), 
{\a\){T,L), and x{T,L), a phase transition likelyoccurs 
for d = 2 and 3, for d = 3 in agreement with Refs-M fTllfl^ . 
However, for d = 1, in spite of the long-range interaction, 
there seems to be no phase transition at finite T. 

This conclusion is confirmed by the behavior of the 
Binder parameter, the fourth-moment ratio, 

Q^{a^f/{a^)^ (7) 
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which is shown in Fig. 4.^ This quantity is directly de- 
rived from the Binder cumulant 1 — ((t^)/(3((t^)^) and 
exhibits similar features: In the ordered phase, where a 
has a bimodal probability distribution, (5 = 1, whereas, 
in case of a normal distribution of u in the disordered 
phase, Q = 1/3. A common intersection point of the 
Q(r) curves for various system sizes indicates a phase 
transition. According to Fig. 4, such points seem to exist 
for d—2 and 3, but not for d~\. 



IV. QUANTITATIVE ANALYSIS BY MEANS 
OF FINITE-SIZE SCALING 

Consider first the one-dimensional case: Careful in- 
spection of Fig. 4 leads to the conclusion that, although 
there seems to be no common intersection point of the 
Q(r) curves, simple scaling QiJ.L) = Q{T/LP) is un- 
likely because of the L dependence of the slope in this 
graph. Thus, to perform a detailed analysis, we numeri- 
cally solved Q{Ta,L) = a for a = 0.45, 0.60, 0.75, and 
0.90. The resulting Ta{L) are presented in a 1/Ta versus 
\nL plot in Fig. 5. 

Figure 5 shows that, with high precision, \/Ta is a 
linear function of InL for all considered A values. (The 
tiny deviation at L = 70 arises from a small systematic 
shift of the Q{T, L) curves for L = An, where n is integer, 
in comparison to the curves with L = 4n + 2. This shift 
vanishes as n — > oo.) Moreover, it is remarkable that the 
slope of the linear function is independent of A within 
the accuracy of the simulations. Finally, with increasing 
L, the L dependence of the maximum temperature Tmax 
of the susceptibility seems to tend to the same behavior 
as TAiL). 

Due to the linear dependences in Fig. 5, Ta{L) and 
Imax(-t') likely vanish as L — > oo. This is confirmed 
by the parallelism of the regression lines, a second argu- 
ment against the intersection of Q{T,L = const) curves 
at finite T. As a consequence of these two observa- 
tions, Q{T,L) is expected to depend only on the com- 
posed quantity z = Tq/T - InL with Tq = 0.2484(8), 
where the uncertainty denotes the 3ct bound of the ran- 
dom deviation. This reduction is confirmed by Fig. 6, 
which shows the corresponding master curve made up of 
Q{T, L = const) curves with L = 40, 100, 280, 700, and 
1400. 

Thus, in the limit L = oo, long-range order should be 
destroyed by thermal excitations at arbitrarily small, but 
finite T although the considered model exhibits a long- 
range interaction. 

One may ask whether the logarithmic relation indi- 
cates that the antiferromagnetic Coulomb interaction is 
a marginal case for d — 1. Thus we performed additional 
simulations for the antiferromagnetic l/r^/^ interaction 
which decays more slowly. These studies yielded simi- 
lar results: Considering samples of up to 1000 sites, we 
could not find a clue of a phase transition at finite T. 
However, on the one hand, 1/Ta rises with increasing L 
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FIG. 5: (Color online) Size dependence of the solution Ta of 
Q{Ta, L) — a for d = 1 in a 1/Ta versus InL plot. From top 
to bottom, the crosses are related to ^ = 0.90 (black), 0.75 
(green / gray), 0.60 (blue / gray), and 0.45 (red / gray). The 
maximum temperatures Tmax of x{T, L) are included as ma- 
genta (gray) circles for comparison. Error bars are omitted 
since they are considerably smaller than the symbol size. The 
straight lines correspond to linear regression in this represen- 
tation for L > 100 in case of Ta, and for L > 700 in case of 
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FIG. 6: (Color online) Check for d = 1 whether Q(T,L) de- 
pends only on the composed quantity z = Tq/T ~ InL where 
To = 0.2484. The symbols A (green / gray), 4- (black), o (red 
/ gray), x (blue / gray), and V (magenta / gray) stand for 
L = 40, 100, 280, 700, and 1400, respectively. 

slightly faster than as a linear function of InL, and, on the 
other hand, small deviations from the parallelism of the 
dependences of 1/Ta on InL for different A are present. 

We turn now to the two- and three-dimensional cases: 
Here, the quantitative evaluation of our simulation data 
consists in a finite-size scaling analysis ^^iSS. However, for 
numerical convenience, we consider the quantities 

(72 = -ln(l-g) (8) 

and 

g3--tan(^(l-1.5g)) (9) 

for d = 2 and 3, respectively, instead of Q. The qd{T) 
have small curvature in the vicinity of the transition. 
This behavior alleviates a precise interpolation. 
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FIG. 7: (Color online) Temperature dependence of g^, related 
to the Binder parameter Q by Eqs. (8) and (9), in the close 
vicinity of the transition for d = 2 and 3. According to in- 
creasing modulus of the slope, the curves refer to L = 16, 24, 
34, 48, 68, 88, and 112 for d = 2, and to L = 8, 10, 12, 14, 
16, 18, 20 and 22 for d = 3. The error bars denote the la 
region. For clarity, data points in the intersection regions are 
omitted. 

Figure 7 shows qd{T) for various L values. For d = 2, 
there clearly is a common intersection point of the curves 
for different L at the critical temperature Tc_2- However, 
for d = 3, only a tendency toward such a behavior is seen 
although the widths of the T intervals for both cases as 
well as the ranges of numbers of sites L'^ are comparable. 

The corresponding corrections to scaling can be taken 
into account to a large extent in a simple way: Follow- 
ing Ref. [33 , we define a size-dependent critical tempera- 
ture Tc^diL) by the demand qd{Tc4{L),L) = qo^d, where 
the QQ d are appropriate constants, which are fixed be- 
low. Scaling of qd{T) curves for different L with respect 
to T — Tc^d{L) yields very good data collapse. Figure 8 
demonstrates this observation showing plots of Q{T,L) 
versus t = ad{L){T — Tc^d{L)) based on the parameter 
values obtained below. In our study, referring to Tc^d{L) 
instead of Tc^di'x) proved to considerably reduce the in- 
fluence of deviations from scaling on the values of the 
critical exponents, which arc numerically obtained from 
samples of finite size. 

Figure 8 simultaneously shows that Q{t) has a substan- 
tial curvature within the crucial region - our finite-size 
scaling analysis is based on the Q intervals [0.75, 0.91] 
and [0.45, 0.80] for d = 2 and 3, respectively. However, 
to reach a high accuracy of the critical exponents, very 
precise ad{L) values are needed so that a broad t range 
has to be taken into account. We approach this nonlin- 
earity problem by considering qd{t) instead of Q on the 
one hand, and by approximating qd{t) by polynomials of 
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-a,{L){T-T,,{L)) 

FIG. 8: (Color online) Scaling check for Q{T,L): The adjust- 
ment of the parameters ad{L) and Tc^d{L) is described in the 
text. Scaling is indicated by the agreement of the nonlinear 
contributions to Q{T, L = const) for various L in these plots. 
For the meaning of the symbols -f, o, and x see caption of 
Fig. 1; A (green / gray) refers to i = 112 and 22 for d — 2 
and 3, respectively. 



third degree, go,d + t + bdt^ + Cdt^, on the other hand. 

We evaluated our qd{T,L) data by a series of regres- 
sion studies of the T and L dependences: First, the only 
weakly i-dependent parameters bd and Cd of the poly- 
nomial ansatz were adjusted, after this the ad{L) and 
Tc,d{L) values were determined. 

Both for the cases d = 2 and 3, the relations ad{L) 
only weakly deviate from power laws, ad{L) oc L^, so 
that they are not graphically depicted here. Based on 
the ad{L) data, the value of the critical exponent of the 
correlation length, v = l/p, can be obtained in two 
ways. First, it can be directly calculated by numerical 
differentiation by means of the midpoint formula utilizing 
h' = (dlnla^l/dlnL)"^. The advantage of this approach 
is a meaningful control of convergence with increasing 
L. Figure 9 shows such results from the consideration 
of pairs of sixth-next and next-nearest neighbors in the 
sequence of sample sizes for d = 2 and 3, respectively. 
Second, v can be determined by means of power-law fits 
taking into account various L intervals. This method 
yields more precise estimates of the exponents. For con- 
sistency, the mean-square deviation of these fits must be 
understandable as resulting from random errors alone. 
Table I presents the most precise results for v, which were 
obtained from the fits safely fulfilling this requirement. 
Additionally, these values are included in Fig. 9. 

This graph implies several conclusions: The ly values 
converge rapidly with increasing sample size so that, for 
d = 3, already the study of samples with L k, 12 yields 
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FIG. 9: (Color online) Approximations of the critical expo- 
nent V versus sample size. The symbols o denote values which 
were obtained by numerical differentiation, see text. The cor- 
responding error bars denote the Icr regions. The thick red 
(gray) lines represent results of power-law fits, where their 
extension visualizes the evaluated L interval. The thin red 
(gray) lines give upper and lower bounds of the correspond- 
ing 3(7 intervals. For comparison, the green (gray) symbols x 
mark the known values of the Ising model with short-range 
interactions for L = oo, see Tab. I. 



TABLE I: Finite-size scaling results for the critical exponents 
of specific heat, mean staggered occupation, susceptibility, 
and correlation length, a, /3, 7, and v, respectively. To re- 
tain numerical precision, we mostly present exponent ratios 
instead of the exponents themselves. The data were obtained 
by two alternative methods, either directly (marked as d) by 
power-law fits (modified by a constant term in the case of 
a/u), or indirectly (marked as i) via Widom and hyperscaling 
relations from power law fits yielding other exponents or expo- 
nent ratios, see text. Values for the Ising model with short- 
range interactional^ are included for comparison. Paren- 
theses and brackets give 3cr random errors and total errors, 
respectively, referring to the last given digit of the value. 



Quantity 


d 


L region 


method 


Coulomb 


s.-r. Ising 


ajv 


2 


28 - 112 


d 


-0.02(4) 


(In) 


a/v 


2 


34 - 112 


i 


-0.03(5) 


(In) 


P/v 


2 


48 - 112 


d 


0.1318(21) 


1/8 


P/v 


2 


48 - 112 


i 


0.129(8) 


1/8 


-y/u 


2 


48 - 112 


d 


1.742(15) 


7/4 


V 


2 


34 - 112 


d 


1.013(25) 


1 


ajv 


3 


10 - 22 


d 


0.09(9) 


0.1 740 [8] 


a/v 


3 


8-22 


i 


0.158(21) 


0.1 740 [8] 


aiv 


3 


14 - 22 


d 


0.506(7) 


0.51820(8] 


P/v 


3 


14 - 22 


i 


0.514(5) 


0.51820(8] 


-y/v 


3 


14 - 22 


d 


1.973(10) 


1.96361(15] 


V 


3 


8-22 


d 


0.633(4) 


0.63012(16] 



0.1296 



'c,3 0.1292 



0.1288 - 
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FIG. 10: (Color online) Relation between a3(I/) and Tc,3{L) 
for qo,3 = —0.200. With increasing as, the points refer to 
L = 8, 10, 12, 14, 16, 18, 20, and 22. The error bar represents 
our extrapolation L —> 00, where the 3a interval is marked. 



v values, which are very close to the results for L w 20. 
Thus the good agreement between our result of a power- 
law fit, 0.633(4), and the known Ising value 0.630, sug- 
gests that the considered model may belong to the Ising 
universality class. This hypothesis is supported by our 
result 0.633(4) clearly differing from the mean-field value 
0.5, from the Heisenberg value 0.71, and from exponents 
of the alternative "nearby" models XY, v = 0.670, and 
self-avoiding walks, v = 0.588, compare Ref. [H. Also for 
d = 2, our v value supports Ising universality and clearly 
excludes mean- field behavior. Finally, according to Fig. 
9, it is not surprising that, for d ~ 3, Overlin et al. got 
the somewhat lower value 0.55 ±0.1 from simulations for 
L = 4, 6, and 8."^ 

In obtaining Tc^d{L) from qd{TcM{L),L) = go,<i, a de- 



viation 6 of qoA from the L 00 limit of the solution of 
qd{T, 2L) = qd{T, L) gives rise to a contribution 6/ad{L) 
to Tc^d{L). We fix qo,d by the demand that this term 
vanishes: go,2 = 1.892(8) and qo,3 = -0.200(19), cor- 
responding to the Q values 0.8492(12) and 0.625(4), re- 
spectively. For d = 2, the critical Q value slightly, but 
significantly deviates from the analytical result for the 
L = 00 limit of the short-range Ising model with peri- 
odic boundary conditions 0.856216."^^ But this is not a 
strong counterargument to the considered model belong- 
ing to the Ising universality class, as it is suggested by 
the value of v: It is known that the critical Q value is a 
very sensitive quantity which depends on the boundary 
conditions, and that "universality of the critical cumu- 
lant holds in a rather restricted sense, when compared 
to universality of critical exponents" Nevertheless, for 
d = 3, our critical Q value perfectly agrees with the value 
for the Ising model with short-range interaction and peri- 
odic boundary conditions, 0.6233(4).'^^ This supports the 
hypothetical Ising criticality infcrcd from the value of v. 

The remaining higher-order corrections in Tc^d{L) orig- 
inate from imperfection of finite-size scaling. Comparing 
several empirical approximations, we observed that, over 
a wide L range, they are almost proportional to ad{L)^'^, 
see Fig. 10. Corresponding extrapolations yield the fol- 
lowing values of Tc,d(oo): 0.103082(9) and 0.128838(17) 
for d — 2 and 3, respectively. The confidence intervals 
include the Scr random errors and cautious estimates for 
the systematic uncertainty of the extrapolation L 00. 

Our Tc.3 value is consistent with the result 0.128±0.005 
in Ref. 7, which was obtained for considerably smaller 
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FIG. 11: (Color online) Approximations ol the critical expo- 
nent ratios 'y/v versus sample size. For the meaning of the 
symbols see caption of Fig. 9. 



p/v 



0.15 
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FIG. 12: (Color online) Approximations of the critical ex- 
ponent ratios P/v versus sample size. The meaning of most 
symbols is defined in the caption of Fig. 9. Here, addition- 
ally, blue (gray) symbols • denote values, which were obtained 
by means of the Widom and hyperscaling relations from 'y/v 
data in Fig. 11, marked by o there. 



samples. This coincidence, together with the approxi- 
mate agreement of the values of v, confirms an assump- 
tion from Sec. II. It shows that, within the accuracy of the 
simulation, our treatment of the surroundings of the sam- 
ple yields the same results as the Ewald summation per- 
formed in Ref. 0- Moreover, comparing with analytical 
theories, we mention that the nonlinear screening theory 
by Pankov and Dobrosavljevic underestimates our "ex- 
act" Tc.s value by 26%)^ whereas a study by Malik and 
Kumar, which utilizes the replica method, overestimates 
Tc.3, by roughly a factor 1.5 when using the random- 
phase approximation, and by more than a factor 3 for 
the Hartree approximationi^^ 

The analysis of x{T,L), {\a\){T,L), and c(r, L) was 
performed similarly to the evaluation of qd{T,L): We 
considered Inx, ln(|cr|), and Inc as functions of t and L. 




FIG. 13: (Color online) Size dependence of the value of (|a|) 
at the critical temperature Tc^d{L) defined in the text. The 
error bars are considerably smaller than the symbol size, thus 
they are omitted. The dashed lines represent the fits given in 
Tab. I. 



For not too large \t\, a,s L ^ oo, scaling implies that each 
of these quantities is decomposable into a sum of two 
functions depending only on t and L, respectively. How- 
ever, for the L regions considered here, this hypothesis 
proved to be well fulfilled only for In^. In the cases of 
Inc and ln(|(T|), there is a clear tendency toward such a 
behavior, but small deviations cannot be neglected. Thus 
we approximated Inx, ln(|cr|), and Inc by polynomials in 
t of third order, taking advantage of universalities in the 
coefficients as far as possible. This regression provides 
precise values for the observables a.t t = 0. Simulta- 
neously, we obtained the confidence intervals taking into 
account the uncertainties in the individual measurements 
of the observables and in the Tc^d{L) values. 

The interpolated xiTcdW^L) and {\a\)iT,^d{L), L) 
were analyzed by means of power-law fits, where propor- 
tionality to L'*'/'^ and L~^/'^ , respectively, was presumed. 
These studies were performed analogously to the deter- 
mination of V. However, while the effective j/i' converges 
quite rapidly with increasing L, see Fig. 11, the determi- 
nation of the limit of /3/^ as i — ^ oo in Fig. 12 is hindered 
by slow convergence: Within a fixed L interval, the rela- 
tive change of the effective f3 is considerably larger than 
the variations of the effective v and 7. The reason of 
this slow convergence is understood by inspection of Fig. 
13. Although this graph shows high-quality power-law 
behavior above L « 30 and L 10 for d = 2 and 3, 
respectively, the uncertainty of the slope is rather large 
since, due to the small relative change of (IctI), small de- 
viations from the power law may considerably shift the 
exponent value. 

Alternatively, the value of (3/v can be obtained from 
the value of 'y/iy utilizing the Widom relation, 2 = 
a + 2(3 + and the hyperscahng relation, 2 — a = di^. 
These equations imply P/u = {d — ^/v)/2. Correspond- 
ing results are included in Fig. 12, as well as in Tab. I, 
additionally to the results of the power-law fits. 
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FIG. 14: (Color online) Approximations of the critical ex- 
ponent ratios a/u versus sample size. The meaning of most 
symbols is defined in the caption of Fig. 9. However, in this 
case, the power-law fit is modified by a background constant, 
see text. Additionally, blue (gray) symbols • denote values, 
which were obtained by means of the hyperscaling relation 
from V data in Fig. 9, marked by o there. The dashed lines 
are included only as guide to the eye in order to demonstrate 
the considerable L dependence of the exponent ratio a/f ob- 
tained by numerical differentiation. 



Table I shows that our j3/p data are very close to the 
exponent ratios for the short-range Ising model, where 
the results of the indirect approach deviate somewhat 
less from the Ising values than the data obtained directly 
by means of power-law fits. 

Compared to the study of {\(j\) and x, the analysis of c 
is more difficult: Exponent values obtained by means of 
numerical differentiation converge only very slowly with 
increasing L, see Fig. 14, and the mean-square devia- 
tions of the power-law fits remain too large. Therefore 
we took into account a background contribution presum- 
ing c(Tc^d{L), L) ^ a + bL^ with p — ajv similar to Ref. 
|39| . Moreover, in order to avoid numerical problems, we 
approximated c(Tc^d(i), by 'a^h{U' — l)/p instead of 
by a + hU" . Results for a/u, which were obtained in this 
way, are included in Tab. I and Fig. 14. 

Alternatively, a/v was calculated by means of the hy- 
perscaling relation utilizing the v results, which were ob- 
tained by numerical differentiation. These a/v values ex- 
hibit a far better convergence than the directly obtained 
a/v data, this is demonstrated Fig. 14. They agree nicely 
not only with the results of the modified power-law fits 
but also with the known values for short-range Ising uni- 
versality, see Fig. 14 and Tab. I. 



der seems to be destroyed already at arbitrarily small but 
finite temperature in spite of the long-range character of 
the interaction. However, for two- and three-dimensional 
systems, continuous phase transitions occur at finite tem- 
peratures. We have determined the critical exponents not 
only for the correlation length, but also for the specific 
heat, for the order parameter staggered occupation, as 
well as for the related susceptibility. 

A survey of the exponent values, which we obtained 
by finite-size scaling, is given in Tab. I. These data have 
to be regarded as effective exponents. Due to the finite- 
ness of L, tiny systematic errors are certainly present, 
presumably the more relevant the smaller the exponent 
value. Unfortunately, our data set is not sufficient for 
a convincing estimate of them. Nevertheless, the expo- 
nents which we obtained directly by means of power-law 
fits (modified by a background constant in case of the spe- 
cific heat) obey the Widom relation, according to which 
the term 2/v — (a/v ^ '2/i/v -I- 7/1^) has to vanish: Pre- 
suming the errors in Tab. I to be independent and ran- 
dom, one obtains -0.01(6) and 0.08(9) for d = 2 and 3, 
respectively, where the errors are given as Scr bounds. Si- 
multaneously, our data satisfy the hyperscaling relation 
which implies the quantity 2/v — a/v — d to equal zero: 
In this case. Tab. I yields -0.01(6) and 0.07(9) for d = 2 
and 3, respectively. Together, Widom and hyperscaling 
relations imply that the expression 2(3 /v -\- ^ / v — d van- 
ishes what can be checked with higher precision since a 
is not included here. For this expression, one obtains the 
values 0.006(16) and -0.015(17) for d = 2 and 3, respec- 
tively, from Tab. I. Thus our set of the critical exponent 
values satisfies all consistency criteria. 

Among the critical exponents, v exhibits the best con- 
vergence with increasing sample size. The values given in 
Tab. I are consistent with the assumption that the stud- 
ied phase transition belongs to short-range Ising univer- 
sality class, in spite of the long-range Coulomb interac- 
tion, both for the cases d = 2 and 3. All other well-known 
universality classes could clearly be ruled out according 
to the values of v. Moreover, the supposed Ising univer- 
sality is supported by the values of the critical exponents 
a, P, and 7 being likewise consistent with the correspond- 
ing critical indices of the Ising model, as well as by the 
critical values of the Binder parameter Q. 

Concluding, in spite of the long-range interaction, the 
Coulomb system described by Eqs. (1) and (2) seems to 
belong to the same universality class as the Ising model 
with short-range interaction. This suggests that screen- 
ing should be highly effective in the ordering process 
and that the lattice Coulomb-glass model might have the 
same critical properties as the random-field short-range 
Ising model. 



V. DISCUSSION 

Summarizing, we have presented a detailed Monte 
Carlo study of the ordering of charges on a half-filled 
hypercubic lattice. For one-dimensional systems, the or- 
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